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ABSTRACT 

We investigate the evolution of the relative angle between the stellar rotation axis and 
the circumstellar disc axis of a star that forms in a stellar cluster from the collapse of 
■ a turbulent molecular cloud. This is an inherently chaotic environment with variable 

accretion, both in terms of rate and the angular momentum of the material, and 
Q . dynamical interactions between stars. We find that the final stellar rotation axis and 

S—{ ' disc spin axis can be strongly misaligned, but this occurs primarily when the disc 

r/i \ is truncated by a dynamical encounter so that the final disc rotation axis depends 

simply on what fell in last. This may lead to planetary systems with orbits that are 
misaligned with the stellar rotation axis, but only if the final disc contains enough 
mass to form planets. We also investigate the time variability of the inner disc spin 
axis, which is likely to determine the direction of a protostellar jet. We find that the jet 
I/"") ' direction varies more strongly for lighter discs, such as those that have been truncated 

lO ■ by dynamical interactions or have suffered a period of rapid accretion. Finally, wc note 

' that variability of the angular momentum of the material accreting by a star implies 

that the internal velocity field of such stars may be more complicated than that of 
. aligned differential rotation. 
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*H i 1 INTRODUCTION of planetary systems in general, because of our special place 

. . . within it, it may also be atypical in some ways (e.g. Beer 

Over the last 20 years or so our picture of the environment , , T , 

■' r et al., 2004). in this paper we consider how an examination 

in which stars form and the way in which they do so has of the degree tQ which steUar rQtation axeg are aligned with 

changed considerably With a few notable exceptions (for their protostellar discs and/or planetary systems can give 

example, Larson 1978), prior to 1990 the process of star for- • c ,. , , . u . , , c m U - • 

* ' " f * us information about the way m which stars form, this is 

mation was treated by forming one star at a time (see, for u ■ n , - A , • . „ , 

■' 6 v ' becoming a relevant consideration, given a recent amount 

example the reviews by Shu, Adams & Lizano, 1987; Mc- r . ■■, c , . ,. , c , 

r ' ' or growing evidence of star-disc misalignment from obser- 

Kee & Ostriker, 2007). The formation of binary stars was , . c ■ , .. • ,. . r , ... , , 

' ' J vations of spm-orbit misalignment of transiting exoplanets 

also considered, but mainly through numerical simulation of (H( s brard et aL; 20 08; Winn et al. 2009a, 2009b; Pont et al., 

the dynamics of rotating rings or bars of material (see re- 2QQQ ^ 2QQ9h . QmQn 2QQQ . Johnson ^ ^ 2QQ9 . N&rita ^ 



al. 2009). 



views by Bodenheimer & Burkert, 200f; Bonnell, 2001). In 
all these scenarios, the symmetries are such that the rota- 
tion axes of any stars formed tend to remain aligned with The fact that many, if not most, stars are members of 
the rotation axis of any surviving circumstellar material or binary or multiple systems has led to the argument that the 
disc. Thus these simple models all predict that protostel- star formation process must be much more dynamical and 
lar discs, and therefore any planetary sytem formed within interactive than the simple pictures described above (e.g. 
them, will be closely aligned with the rotation axis of the Pringle, f989, 1991; Clarke & Pringle, 1991). Early multi- 
central star. This expectation is reinforced by observation of plicity surveys of pre-main-sequence stars in nearby star- 
the one system for which we have good data - the Solar Sys- forming regions confirmed that most low-mass stars in the 
tem (Tremaine, 1991). We need to be aware, however, that solar neighbourhood form in multiple systems (e.g. Leinert 
while the solar system may be a good guide to the formation et al., 1993; Ghez, Neugebauer & Matthews 1993; Simon et 
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al., 1995). Indeed most stars within the Galaxy appear to 
form in clustered environments (see the reviews by Lada & 
Lada, 1995, Clarke, Bonnell & Hillenbrand, 2000) and there 
is evidence that disc sizes and presumably disc lifetimes are 
reduced in binary as compared to single stars (Cieza et al., 
2009). The theoretical expectation is that close encounters 
in the early stages of formation when stars have substantial 
discs lead to substantial disc truncation and moreover, that 
the tidal effects of such encounters can also lead to the outer 
parts of the disc becoming strongly misaligned or warped 
(Heller, 1995; Hall, Clarke & Pringle, 1996; Larwood, 1997; 
Pfalzner, 2004; Moeckel & Bally, 2006). Such interactions 
imply that there is a strong possibility that the stellar rota- 
tion axis (determined by the angular momentum of material 
accreted early which forms the stars) is not aligned with the 
disc rotation axis (which might correspond to the angular 
momentum of material accreted later which forms the disc, 
and eventually planets). 

Recent numerical simulations of the star formation pro- 
cess (Bate, Bonnell & Bromm, 2002a,b, 2003; Bate & Bon- 
nell 2005; Bate 2005, 2009a,b) have considered the formation 
of stars within clusters, starting with initial turbulent con- 
ditions for the gas which are designed to mimic those seen 
within star forming molecular clouds. The picture which re- 
sults involves turbulent and chaotic motions of both gas and 
stars, with disc fragmentation, competitive accretion and 
close dynamical interactions all playing a role. Close star- 
disc encounters occur throughout the star formation process 
at a variety of orientations. From this perspective it might 
seem reasonable to conclude the star-disc misaligment would 
be the norm rather than the exception. Here we investigate 
this possibility. 

Because of the large range of scales involved it is not 
feasible to follow the evolution of gas starting from the tur- 
bulent molecular cloud (tens of parsecs) down the the sizes 
of protostellar discs (tens of au) and further to the accre- 
tion of matter onto the protostars themselves (a few solar 
radii). Numerical simulations of the star formation process 
begin with a turbulent molecular cloud but truncate the 
computation typically at a resolution scale of several au by 
means of sink particles. In this paper we are interested in 
what happens to the gas when it falls within the sink par- 
ticles. We are therefore constrained to carry out a highly 
idealised computation using a strongly simplified, and not 
necessarily self-consistent, set of assumptions. We consider 
the evolution gas falling into a sink particle taken from a 
particular simulation. In Section 2 we describe the simula- 
tion. In Section 3 we describe the simplifications we have 
made in order to model the evolution of the gas when it has 
fallen into the sink particle from the formation of the sink 
particle (« 10 -3 Mq) until after accretion onto the sink par- 
ticle has finished and the 'star' is more or less fully formed 
with mass « 0.2 Mq. In Section 4 we present the results of 
our calculations. We discuss the implications in Section 5. 



2 SIMULATION DATA 

We are interested here in illustrating the growth of a partic- 
ular protostellar core. The data we use comes from a com- 
putation by Bate, Bonnell & Bromm (2003) which describes 
the formation of a cluster of stars from a turbulent self- 



gravitating parcel of dense (molecular) gas. The details of 
the numerical procedures are described fully in that paper, 
and here we draw attention to those details of particular 
relevance to the calculations below. 

The calculation was performed using a 3D smoothed 
particle hydrodynamics (SPH) code. The smoothing lengths 
are variable in time and space, subject to the constraint that 
the number of neighbours for each particle must remain ap- 
proximately constant at iV no igh = 50. The initial cloud is 
spherical and uniform in density with a mass of 50 Mq and 
diameter of 0.375 pc (77,400 au). At the initial tempera- 
ture of 10 K, the thermal Jeans mass is 1 Mq. The free 
fall timescale of the cloud is ts = 1-90 x 10° yr. An initial 
suspersonic turbulent velocity field is imposed on the cloud 
with rms Mach number M = 6.4. The field is divergence-free 
and random Gaussian with a power spectrum P(k) oc fc~ 4 , 
where k is the wavenumber. The local Jeans mass is resolved 
throughout the calculation. The number of particles used is 
3.5 x 10 6 , which means that each particle has a mass of 
1.43 x 10 -5 M Q . 

The calculation uses a barotropic equation of state that 
is isothermal up to densities of 10 -13 g cm~ 3 , but with the 
temperature of the gas increasing at higher densities. In this 
way, the calculations mimic the opacity limit for fragmen- 
tation (Low & Lynden-Bell 1976; Rees 1976). Despite this, 
because of computational constraints, it is necessary to ex- 
cise the high density protostellar cores from the bulk of the 
computation. This is done by inserting a sink particle (Bate 
et al., 1995) whenever the central density of a pressure- 
supported fragment exceeds a certain density, p cr it = 10 -11 
g cm~ 3 . The sink particle is formed by replacing the SPH 
gas particles contained within a radius i? S mk = 5 au of the 
densest gas particle by a point mass with the same mass 
and momentum. Any gas that later falls within this radius 
is accreted by the point mass if both (a) it is bound and 
(b) its specific angular momentum is less than that required 
to form a circular orbit at radius i? s ink from the sink parti- 
cle. Sink particles interact with the gas only via gravity and 
by accretion. Sink particles interact with each other only 
by gravity. The gravitational acceleration between two sink 
particles is softened within a distance of 4 au. Thus sink 
particles cannot merge with each other. The typical initial 
mass of a sink particle in this calculation is w 10~ 3 Mq. 



3 EVOLUTION EQUATIONS 

As we have seen, because of resolution constraints, in the 
numerical simulations of chaotic star formation the 'stars' 
formed are represented by sink particles. Because of this, 
from the simulations we only have very basic information 
about material being accreted onto the protostellar cores. 
We have information about the rates of accretion of both 
mass and angular momentum onto the sink particles, but 
note that the accretion is discretized into that of individual 
SPH particle masses of « 10 Mq. The sink particle radius 
is set at 5 au., which is slightly less than the peak in the 
distribution of binary star separations for low mass stars. 
We also have information about the distance to the nearest 
neighbouring sink particle as a function of time. It is usually 
the case that for some of the time this distance is less than 
the radius of the sink particles. 
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Since the information about the accretion process onto 
sink particles is so limited, we are strongly constrained in 
how we are able to model what might happen to the gas 
when it is accreted within the sink particle radius of 5 au. We 
proceed by making the simplest assumptions we can about 
the gas flow within the sink particle, noting that some of 
the implications of these are of necessity at times not self- 
consistent. 

Within the sink particle, we model the gas flow as a 
twisted accretion disc around a central mass. The evolu- 
tion of such a disc is computed in the manner described by 
Pringle (1992). The disc surface density is given as a func- 
tion of time and radius by E(r, i). Each disc annulus, at 
radius R, is assumed to rotate with angular velocity ii(r, t) 
around the central mass, and to have a spin in the direction 
of the unit vector 1(_R, t). Thus the disc locally has an angu- 
lar momentum density given by L(7?,i) = 'SR 2 Ql. We shall 
also assume (cf. Lin & Pringle 1990) that any matter being 
added to the disc is added at the radius corresponding to its 
specific angular momentum, so that it arrives already on a 
circular orbit and in centrifugal balance. 

Then (Pringle 1992) the evolution of the angular mo- 
mentum density is given by 



dt 
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(1) 



Here Q,' — dQ/dR, v\ is the shear viscosity normally associ- 
ated with accretion discs and v 2 is associated with straight- 
ening out the out-of- plane motions (i.e. the warp or twist). 
The source term S(R,t) allows for the addition of material. 

When the sink particle forms it has a mass, typically 
around ~ 10 -3 Mq, which is much less than the eventual 
mass of the star. Thus during the accretion process we ex- 
pect the overall mass of the star+disc system to grow by a 
factor of around 10 2 — 10 3 . To calculate the angular velocity 
Q(R) we need to take account of the central time-varying 
force contributions from both the central star, mass M*(i) 
and the disc itself. Although an exact determination of Q(R) 
could be made (e.g. Bertin & Lodato 1999), in general we 
shall find that the disc mass is much less than of the central 
star. Thus for computational convenience (cf. Lin & Pringle 
1990) we adopt the approximation 
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where 



M(R,t) = M„(t) + S(i?, t) 2nR dR 
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Jo 



(2) 



(3) 



In order to estimate the magnitudes of the viscosity 
terms it is strictly necessary to compute the disc structure 
at each radius. For protostellar discs this is not a straightfor- 
ward matter and the physical processes relevant at various 
radii are still a matter for debate (see for example Lodato 
& Rice 2004, 2005; Kratter, Matzner & Krumholz, 2008; 
Terquem 2008). What would be required first is a detailed 



time- dependent model of the evolution of the disc surface 
density distribution, similar to that carried out by Lin & 
Pringle (1990), but taking more recent concepts and un- 
derstanding into account (see, for example, Armitage, Livio 
& Pringle 2001; Rice & Armitage 2009; Zhu, Hartmann & 
Gammie 2009). In addition, a secondary, but nevertheless 
important, consideration would be to take account of how 
such discs might respond to misalignments - for example, 
the low viscosity regions (the dead zones) might transfer 
warp in a wave-like, rather than in a diffusive, manner, 
and gravitational bending torques might play a role in disc 
regions where self-gravity is important (Papaloizou & Lin 
1995). 

With these issues in mind, and in order to simplify 
matters in this initial investigation, we make the follow- 
ing assumptions. We assume that for each viscosity, i.e. for 
i = 1,2, 



(4) 



Here H is the disc semi-thickness, and we shall assume for 
convenience that H/R = 0.1, independent of radius (cf. Bell 
et al., 1997; Terquem, 2008). Thus Qi is the usual Shakura 
& Sunyaev (1973) viscosity parameter and we take typically 
(cf. King, Pringle & Livio, 2007) 



0.02, 



(5) 



although we also consider one case with qi = 0.2. Typically 
Q2 ^> c*i (Papaloizou & Pringle 1983) and we shall take 
(Lodato & Pringle 2007) 



Q2 = 2. 



(6) 



As far as the evolution of the system is concerned, what 
matters most are the relative sizes of the various physical 
timescales. Thus to zeroth order, the details of the viscous 
processes are less important than the magnitude. This gives 
rise to the viscous timescale for mass flow through the disc 
which is then given by 



R 



3/2 



.5 J 



The timescale on which warp is propagated through the disc 
is less than this and is 



90 



V 2 ) V0.LR 



Mr, 



m y 1/2 ^- r y /2 



5 au/ 



yr- 



(8) 



3.1 Numerical details 

We take a fixed radial grid with inner radius R ln = WRq. 
The choice of an outer radius of the grid presents us with a 
consistency problem. Some of the sink particles in the calcu- 
lation have resolved discs outside of the sink particle radius. 
The interactions between these discs and the discs which we 
model within the sink particles is not accounted for. This 
is problematic because the timescale governing the propa- 
gation of warps is much shorter than the viscous evolution 
timescale. Therefore, the direction of the spin axis in the disc 
inner parts is affected by that of the spin axis in the disc 
outer parts, which may contain most of the mass. However, 
even if the simulation of Bate et al. were to be performed 
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again, there is no way to provide a continuous and consis- 
tent link between what is inside and what is outside the sink 
particle. When sink particles were originally invented, even 
the construction of boundary conditions that tried to min- 
imise the unphysical effects on pressure and viscosity that 
were introduced outside the sink particle accretion radius 
proved extremely difficult to implement (Bate 1995; Bate, 
Bonnell & Price 1995) . There is certainly no way to link the 
propagation of warps across the accretion radius. Instead, 
for the calculations presented here, we have chosen one of 
the many sink particles that does not have a resolved disc 
outside of the sink particle radius. 

However, even for those sink particles not surrounded 
by resolved discs, much of the material accreted onto the 
sink particle has an angular momentum which would put 
it on a circular orbit at just inside the sink particle radius 
of 5 au. Thus we expect the disc evolution to lead to disc 
expansion outside 5 au. This is of course not permitted in 
the original calculation. However, it would be unrealistic to 
choose the outer disc radius at 5 au and so to simply re- 
move all disc mass which expands beyond 5 au as that leads 
to a significant fraction of the material accreted onto the 
sink particle not being accreted by the central star. Instead 
it would be lost back to the computational domain. Thus 
whatever choice we make leads to an inconsistency. As a 
compromise we choose the outer disc radius to be R ou t = 50 
au and note that we therefore allow the disc to expand out- 
side the radius of the sink particle i? s i n k = 5 au. 

We employ an explicit first-order finite difference 
scheme, as detailed in Pringle (1992), to solve the time- 
evolution of equation [T] We take 120 logarithmically spaced 
grid points, Thus the ratio of radii of neighbouring grid 
points is 1.06 so that dR/R = 0.06. The timestep is ad- 
justed to satisfy the usual numerical stability criteria. For 
most of the time it is set by the diffusion time-scale across 
the innermost radial zone. 

At each time step the following procedure is carried out. 

1. Mass is added to the disc at the appropriate radius. 
This is the radius at which the specific angular momentum of 
the added matter equals that of the local disc material. The 
largest radius at which matter can be added corresponds to 
the the radius of the sink particle 7? s ink = 5 au which is a 
factor of 10 smaller than the outer disc radius at all times. If 
the added material has sufficiently small angular momentum 
that it lies within the inner disc radius Ri n ~ WRq then 
the mass and angular momentum is added to the central 
star. Because the mass of individual SPH particles are quite 
massive compared to the mass of the disc we spread the 
accretion of mass onto the disc in time by assuming that 
at each time t the accretion rate is a constant given by the 
mass of an SPH particle divided by the time between its 
arrival and that of its predecessor. 

2. At this stage because mass has been added, and be- 
cause disc evolution has taken place (see step 3), the angular 
velocity profile of the disc must be adjusted. First the an- 
gular velocity profile is updated using equation [2] Then the 
amount of material in each grid zone is adjusted to restore 
centrifugal balance. To achieve this we move material in- 
wards in the manner described by Bath & Pringle (1981; 
see also Lin & Pringle, 1990). This process conserves an- 
gular momentum to machine accuracy, but conserves mass 
only to the order of the numerical scheme. 



3. The evolution of the disc angular momentum distri- 
bution is computed from equation [T] The numerical scheme 
is written in such a manner that it conserves angular mo- 
mentum to machine accuracy, but conserves mass only to 
the order of the scheme. At the inner edge of the grid we 
apply a zero-torque boundary condition. Mass and angular 
momentum are freely accreted onto the star, and the stellar 
mass M* (t) and angular momentum L* (t) are appropriately 
updated. At the outer edge of the grid we allow angular mo- 
mentum (and mass) to flow freely outwards. 

3.2 Mass conservation 

We have noted that the although the various steps in the 
numerical scheme conserve angular momentum to machine 
accuracy (except for the small loss at the outer disc edge), 
they do not conserve mass. We show in the Appendix that 
for Step 2 this leads at each timestep to a systematic over- 
estimate of the mass added by a factor of order (dR/R) 2 . 
Similar considerations apply to Step 3. 

3.3 Interaction with other sink particles 

At each time we have the distance to the nearest sink par- 
ticle -Rnext(i) and have noted that it is often the case the 
-Rnext < Rsink- We are only able to take the interaction with 
the neighbouring sink particle into account in an idealised 
fashion. To do so we take note of the finding (e.g. Hall et al, 
1996) that the main effect of the interaction between a disc 
and the fly-by of a point mass is to unbind the outer disc re- 
gions and so to truncate the disc. Here we shall mimic this 
interaction by removing all disc material which is at radii 
R < 0.5-Rncxt- We also note that an encounter is likely to 
perturb the angular momentum of the disc, particularly near 
the truncation radius. We are forced to neglect any such ef- 
fect, but we note that a typical dynamical encounter is also 
associated with the accretion of new material into the sink 
particle and that the addition of this new material, whose 
angular momentum is likely to be uncorrelated with that 
of the existing star and disc, is likely to be more important 
than the tidal effect of the encounter on the pre-existing disc. 
Indeed, as will be seen in the calculations presented below, 
when truncation of the disc through interactions with other 
sink particles is included, the stellar and disc axes ending up 
strongly misaligned because the angular momentum of the 
new material is unrelated to the star's angular momentum. 
If anything, including the perturbation suffered by the pre- 
existing disc during an encounter would typically serve only 
to increase the degree of misalignment, not decrease it, and 
therefore would not alter the fundamental result. 

3.4 Spin directions 

Since we know the total angular momentum accreted onto 
the central star we can compute the direction of that vec- 
tor in space. We use an arbitrary fixed coordinate system 
and denote the space direction of the stellar spin vector in 
spherical polar coordinates as (0*(t), <j)*(t)). If and when the 
star is able to smooth out any internal differential rotation 
this would represent the stellar rotation axis. 

We also measure the spin direction of the disc. We take 
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this to be represented by the spin direction (9 d (t), <f>d(t)) of 
the annulus close to the inner disc edge at » 3.2i?i„ = 32i? . 
We choose this radius as it gives a reasonable estimate of 
the direction in which any jet might be launched, since the 
velocity of such jets indicates that they must be launched 
from close to the central star (e.g. Pringle, 1993; Livio, 1997; 
Price, Pringle & King 2003) . In addition, by the end of the 
calculation at a time of 5 x 10 4 yr, taken to be well after 
the last accretion event, the disc is essentially coplanar, and 
so this direction represents the plane of the protostellar disc 
and any subsequent planetary system. 



4 RESULTS 

In this section we consider the accretion history of one par- 
ticular sink particle taken from the 50 which formed during 
the calculation. Because of the inherent uncertainties and 
inconsistencies of dealing with closely interacting sink par- 
ticles we choose to follow a sink particle which ends up as a 
single star. The total mass accreted onto the sink particle is 
shown as a function of time from the formation of the sink 
particle in Figure [T] (black solid line). It can be seen that 
accretion continues until a time of 8 x 10 11 s = 2.5 x 10 4 
yr. The total mass accreted is 4 x 10 32 g « 0.2 M . Also in 
Figure [1] we show the distance to the nearest sink particle 
Rnext as a function of time (red short-dashed line) . The sink 
particle is the second sink particle that formed in the orig- 
inal hydrodynamical calculation and starts its evolution as 
a companion to the first sink particle formed in the calcula- 
tion. As can be seen from Figure [1] for some 2 x 10 4 yr these 
two sink particles form a binary with separation of around 
20 au. After that, the binary encounters a tighter binary 
consisting of the 3rd and 10th objects to form in the cal- 
culation. After a chaotic interplay lasting about 3000 years, 
during which the distance of closest approach occasionally 
become as small as 1 au, the sink particle we are considering 
is ejected from the system (at the time of 2.5 x 10 4 yr in Fig- 
ure Q]) leaving objects 3 and 10 as a very tight 1 au binary 
(that survives until the end of the calculation) with the first 
sink particle as a wide companion to the tight binary with a 
separation of about 60 AU. This wider companion does not 
survive until the end of the calculation - it is unbound dur- 
ing another dynamical encounter 5000 years later. Accretion 
onto the sink particle that we consider in this paper comes 
to an abrupt end when it is ejected from the multiple system 
and is left as a single star. The computation continues until 
a time 5 x 10 4 yr has elapsed. 

As can be seen from Figure Q] the accretion rate is not a 
constant. Moreover, the accretion onto the sink particle oc- 
curs from a variety of directions. We have seen that it is not 
possible to make a set of fully consistent assumptions about 
how to treat the material that falls into the sink particle. For 
this reason we shall consider two extreme possibilities. First, 
we follow the evolution of the disc ignoring any interaction 
with neighbouring sink particles. In this way we model a 
star forming in the centre of a protostellar disc which is 
being fed material from a variety of directions. Second, we 
take account of the interaction of neighbouring sink parti- 
cles which pass through the disc by simply truncating the 
disc in the manner described in Section 13.31 and removing 
the truncated matter. This is, of course, not fully realistic 
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Figure 1. Total mass accreted by the sink particle (solid black 
line), distance to the nearest other star in au (short-dashed red 
line) , and mass accreted by the central star (long-dashed blue line) 
as functions of time. Although the nearest other star approaches 
within a few au, the effect of the encounters are ignored for the 
evolution of the disc in this case. 
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Figure 2. The mass accretion rate of the central star (solid black 
line) and the orientation angles of the stellar rotation axis and the 
rotation axis of the inner disc as functions of time. The stellar 
rotation axis angles (#» = [0, 180], <f>* — [—180,180] degrees) are 
given by short-dashed lines in red and blue, respectively. The 
inner disc axis angles (S^,ipd) an d given by long-dashed lines. 
After the first 5000 years or so the stellar rotation and inner disc 
axes are well aligned. In this case, the effect of encounters with 
other stars on the disc evolution is ignored. 
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Figure 3. Same as Figure [T] except that the disc viscosity is an 
order of magnitude larger at a± = 0.2. Note that the mass passes 
through the disc more quickly than in the low viscosity case and so 
the central star grows in mass more quickly and because the disc 
drains quickly, accretion onto the central star stops soon after the 
accretion into the sink particle ends. The disc does not maintain 
a significant reservoir of material. 
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Figure 4. Same as Figure [2] except that the disc viscosity is an 
order of magnitude larger at cti = 0.2. Note that the mass accre- 
tion rate on to the central star is higher than the low viscosity 
case earlier, but the accretion onto the central star dies off soon 
after the mass accretion into the sink particle ceases. In addition, 
the stellar and inner disc axes are more misaligned than in the low 
viscosity case. Because accreted material passes quickly through 
the disc the mass of the disc is low and its angular momentum 
depends more sensitively on the material most recently accreted. 



as in reality the material should be put back into the SPH 
computational domain. There it would interact with other 
circumstellar material and would stand a chance of being 
re-accreted. 



4.1 Without sink particle interaction 

We first consider the evolution of the system without tak- 
ing account of any interaction with the neighbouring sink 
particle. As we mentioned above, although not strictly self- 
consistent, this serves to illustrate the general behaviour of 
a star-disc system which is accreting material with a wide 
range of angular momenta. 

We consider the case when the disc viscosity is given 
by Qi = 0.02. In Figure [2] we show the accretion rate onto 
the central star as a function of time (solid black line) and 
the time-evolution of the directions of the stellar rotation 
axis(#*,0*; short-dashed lines) and of the rotation axis of 
the inner disc (6d,(f>d; long-dashed lines). At the end of the 
computation the stellar mass is A/* = 0.113 Mq and the 
disc mass is M d = 3.4 x 10" 2 M = 0.3 M» (self-gravity is 
likely to play a role in providing the effective viscosity of the 
disc in this particular case; Lodato & Rice 2004, 2005). Thus 
about a quarter of the mass accreted onto the sink particle 
has been lost at the outer disc edge. This is to be expected 
since much of the matter is accreted with an angular mo- 
mentum corresponding to a centrifugal radius ~ -R S mk and 
angular momentum conservation leads to the disc expand- 
ing beyond this radius. The fraction of matter required to 
carry away the angular momentum to radius Rout is approx- 
imately l/R s ink/-Rout ~ 0.32. 

For Qi = 0.02 the accretion timescale (Equation]?! from 
where most the matter is accreted into the disc (< 10 4 yr 
at radii R < ii s i n k = 5 au) is comparable to the typical 
timescale for the growth of stellar mass (» 2xl0 4 yr). Thus, 
there is still a non-negligible amount of mass in the disc 
at the end of the computation. At that time the accretion 
rate onto the star is M « 3.8 x 10 -7 Mq yr -1 . Comparing 
these values with the conventional picture of the evolution of 
low-mass stars, our modelling corresponds to the Class 0/1 
phases of protostellar evolution and ends at a stage compa- 
rable to the beginning of the Class II or T-Tauri phase, both 
in terms of disc mass and stellar accretion rates. During the 
evolution, the stellar rotation axis and the disc spin axis 
both converge quite quickly to their final values. Because 
the disc evolution timescale is comparable to the stellar ac- 
cretion timescale, although its spin direction is somewhat in- 
fluenced by late accretion of material (especially when that 
material is accreted at radii close to the inner disc radius 
where the disc spin is measured) it does not differ substan- 
tially from that of the star. At the end of the computation 
the angular distance between these two axes is 4.2°. 

In Figures [3] and [4] we consider an identical case, except 
that the disc viscosity is increased by an order of magnitude 
to ct\ = 0.2. In this case the accretion timescale is much 
less than the timescale for the growth of the stellar mass. 
Thus, the gas passes through the disc quickly, the central 
star grows more quickly in mass (Figure [3]), the disc mass 
is lower at any given time, and significant accretion from 
the disc onto the central star ceases at essentially the same 
time as the accretion into the sink particle ends (Figure 
[4]). Because the disc mass is lower, the spin of the inner 
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Figure 5. Same as Figure[T] except that the disc is truncated by 
the encounters with the nearest other star. Note that the mass 
accreted by the central star is dramatically reduced when the 
effects of the encounters are taken into account because the disc 
is repeatedly truncated. 
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disc is perturbed more strongly by incoming material and 
so the stellar and disc spins are generally less aligned with 
each other than in the low- viscosity case. At the end of the 
computation, the stellar mass is M, = 0.122 Mq, the disc 
mass is only M d = 3.0 x 10~ 5 Mq = 2.5 x 10~ 4 M*, the 
accretion rate on to the star is only M » 2.3 x 10 -9 Mq 
yr -1 , and the angle between the stellar and disc rotation 
axes is 21.3°. 



4.2 With sink particle interaction 

We have seen in Figure [T] that thoughout the time accre- 
tion is taking place there is a nearby sink particle in close 
attendance. For the first 2 x 10 4 yr or so the sink parti- 
cle is sufficiently far away (i? nex t ~ 20 au) that it does not 
interfere strongly with the evolution of the disc, although 
its occasional sorties to within 10 au, imply (under our as- 
sumptions of disc truncation) that loss of angular momen- 
tum from the outer disc regions can be significant. For the 
final « 5000 yr of accretion, the neighbouring sink parti- 
cle regularly approaches within 4 au, before becoming un- 
bound. Thus during late accretion the disc is severely trun- 
cated. This implies that the material which is in the disc at 
the end of the computation will have arrived only shortly 
before the neighbouring sink particle vanished and accre- 
tion terminated. For comparison with the above we take 
Qi = 0.02. In Figure [5] we show the amount of material 
which has accreted onto the central star as a function of 
time (long-dashed blue line). At the end of the computation 
the stellar mass is M* = 1.1 x 10 -2 Mq and the disc mass is 
M d = 2.2 x 10~ 4 M = 2.1 x 10" 3 M,. The accretion rate is 
M « 10 -10 Mq yr -1 . Thus the effect of a nearby companion 
has been to severely decrease the amount of mass accreted 
onto the central star. This is, of course, a direct result of our 
somewhat extreme assumptions about truncation. 

In Figure [6] we show the angular coordinates of the stel- 
lar angular momentum vector (#», 0*) and of the inner disc 
spin vector (#d,0d) as functions of time. In this case it is 
apparent that while the companion sink particle keeps its 
distance, the spin axes evolve more or less in tandem as be- 
fore, except that since the disc is now less massive it is more 
easily buffeted by incoming material. Thus we see that the 
disc spin axis can be quite variable throughout the accre- 
tion process (The points here are separated by 10 year in- 
tervals.). Moreover, during the brief final episode of chaotic 
close encounters disc spin axis becomes erratic as the disc is 
stripped. The upshot of this is that the final disc spin axis 
is unrelated to that of the star, with the angular distance 
between them being 122°. 



Figure 6. Same as Figure[2] except that the disc is truncated by 
the encounters with the nearest other star. Note that the mass 
accretion rate on to the central star is dramatically reduced when 
the effects of the encounters are taken into account because the 
disc is repeatedly truncated. In addition, throughout most of the 
evolution the stellar and inner disc axes are significantly mis- 
aligned. When the star is finally ejected as a single object this 
misalignment persists. 



5 DISCUSSION 

We have considered in an approximate and idealised man- 
ner what happens to gas when it falls into a sink particle 
in a numerical simulation of a cluster of stars. Material is 
accreted onto the sink particles at a variable rate and from 
a variety of directions. The aim of our calculation has been 
to try to understand how the direction of the stellar rota- 
tion axis might evolve, how the direction of the spin axis of 
the inner disc (presumably related to the direction of any 
accretion-induced jet) might evolve, and to what extent, if 
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any, the final stellar spin direction and the final disc plane 
(and hence any resulting planetary orbits) might be related. 

We have limited ourselves to considering the evolution 
of a single sink particle from a particular simulation. In addi- 
tion because of the limitations of the data available, and the 
nature of assumed accretion by sink particles, the nature of 
the computations we are able to carry out are of necessity 
not fully self-consistent and contain various idealisations. 
Thus our efforts here should be regarded as an illustrative 
exploration of the kind of things that might happen rather 
than definitive predictions. Nevertheless we feel that are able 
to draw some general conclusions which indicate that further 
consideration of this problem will be worthwhile. 

The evolution of the direction of the stellar spin axis 
is probably reasonably well modelled. It illustrates how the 
spin axis might evolve given that the angular momentum 
of the material from which the star is forming is signifi- 
cantly variable in direction. The total stellar spin direction 
converges quite quickly to its final value. This is because 
the inner disc direction is controlled quite closely by the 
plane of the outer disc (equation [SJ and thus the average 
spin direction of material arriving in the disc is communi- 
cated efficiently to the star. Thus to a first approximation 
the stellar spin direction just reflects the sum of everything 
that arrives and once most of the mass has been accreted it 
is essentially fixed. What we do not know, however, is how 
the direction of the surface spin of the star relates to its 
total angular momentum. Since these stars are fully connec- 
tive for much of the time it may be that the surface follows 
the total quite well. But we have here a significant difference 
from the usual considerations in that the accretion does not 
all occur axi-symmetrically and thus the internal mixing of 
angular momentum in the star is much more complicated 
than what is usually considered. This has implications for 
dynamo driving - the driving is likely to be quite severe as 
the star sorts out the fact that different shells started ro- 
tating about different axes. Thus the usual relationships be- 
tween surface rotation and magnetic activity might need to 
be modified. Thus for very young stars we might expect that 
the dynamo activity depends as much on accretion history 
as on its surface rotation rate. Further, as stars with masses 
above 0.4 Mq approach the main sequence, their cores be- 
come radiative and so able to decouple from the convective 
outer layers. It might be that these central radiative regions 
can retain some memory of the variable spin direction of the 
formation process. 

As a measure of the direction of disc spin we have 
taken the direction of the angular momentum vector of 
the innermost disc radii. We do this for two reasons. First, 
while accretion is active, this is likely the direction of any 
accretion-driven jet. And, second, because the warp smooth- 
ing timescale (equation [SJ is relatively short, once accretion 
onto the disc has ceased this is a good measure of the spin of 
the final disc as a whole. Given the limitations of the analy- 
sis discussed above, we limit ourselves to general comments 
about our expectations for the evolution of the direction of 
disc spin. We take as two possibly extreme examples of the 
nature of the accretion process the case when the disc retains 
a substantial fraction of the total mass (Figures [T] [5} and 
the case when the disc contains little mass, either because 
it is severely truncated by interactions with a neighbouring 



star (Figures [5j [6} or has undergone rapid accretion (Figures 

HE}. 

Since the accretion process is variable both with regard 
to mass flux and with regard to vectorial angular momentum 
direction, the disc spin can in both cases be variable during 
the process of star formation. While the angular momentum 
of the accreted material is relatively high, so that the mass 
is added at relatively large radii, the accretion process is 
smoothed by the diffusive effects of the disc and the direction 
changes fairly smoothly. However, as is particularly evident 
in the truncated case for which the disc mass is low, there 
are times when the angular momentum of the material being 
accreted is sufficiently low that material is added at small 
radii and the disc direction changes rapidly. Thus during 
the period of close interaction when the disc is continuously 
truncated by repeated close encounters, the disc mass is low 
and so the smoothing effect of the disc is reduced. At these 
times the disc spin direction can vary quite strongly and 
quite rapidly. This also occurs without truncation if the disc 
mass is low because the viscosity is high. It has long been 
recognised that the outflows from young stars - the Herbig- 
Haro flows and jets - can provide probes of early stellar 
evolution (Reipurth & Bally, 2001). It was argued early on 
(Stahler 1994) that the broad molecular flows seen around 
young objects are most likely driven by entrainment by high 
speed jets emanating close to the central star and it was 
also clear that a jet with variable direction is more able to 
entrain material efficiently (Stone & Norman, 1994; Smith 
et al., 1997). In addition there is considerable evidence that 
jets and collimated outflows vary both in outflow rate and 
direction (see for example Cunningham, Moeckel & Bally, 
2009; and the reviews by Bally, Reipurth & Davis, 2007; 
Bally, 2007). In addition we have seen that accretion of mass 
and angular momentum becomes particularly erratic during 
close interactions with companions, in line with the ideas of 
Reipurth (2000). 

We have seen that once accretion onto the disc has 
ceased, the disc settles down to its final planar configura- 
tion well before the disc is able to drain onto the star. This 
occurs because typically U2 S> V\. The correspondence be- 
tween the spin axis of the final disc and the spin axis of 
the star depends on the extent to which the final disc has 
memory of the direction vector of the total anglar momen- 
tum accreted onto the disc, and thence communicated to 
the star. Thus if the disc remains essentially undisturbed 
throughout accretion process, the stellar and disc rotation 
axes are reasonably well aligned. However, if the outer parts 
of the disc (where most of the angular momentum resides) 
are lost, then the final disc direction corresponds to the spin 
axis of the most recently accreted material. This axis might 
or might not be related to the spin of the star. Similarly, 
the action of removing the disc itself (i.e. a dynamical en- 
counter with another object) may misalign the remnant disc 
from the stellar spin axis. 

It is interesting to speculate about the possible depen- 
dence of misalignment on stellar mass. For example, in the 
hydrodynamical calculation presented by Bate et al. (2002, 
2003) and subsequent calculations brown dwarfs are pro- 
duced when protostars that have recently begun to form 
suffer dynamical encounters and are kicked out of the cloud 
terminating their accretion before they have been able to 
accrete to stellar masses. In this case, the brown dwarf may 
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be left with a misaligned disc if either the disc and stellar 
rotation axes had not had time to converge (because the ob- 
ject had recently formed) or if the dynamical encounter itself 
changed the discs rotation axis. Because the object has been 
ejected from the cloud it is unlikely to subsequently accrete 
material with different angular momentum. Conversely, in- 
termediate and massive stars in such hydrodynamical cal- 
culations typically suffer many dynamical encounters with 
other objects (e.g. Bate et al. 2003; Bate & Bonnell 2005; 
Bonncll & Bate 2005) . These encounters often leave the star 
embedded in the dense gas and, thus, repeatedly destroy 
their discs but also allow them to reform their discs from 
freshly accreted material. Again, this may lead to discs be- 
ing misaligned with the stellar rotation axis if the majority 
of the stellar mass was accreted before the final encounter 
and growth of a new disc. 

From the observations of transits of hot Jupiters, re- 
cent measurements of the angle between the planetary or- 
bital plane and the stellar rotation direction show that in a 
substantial fraction of cases the orbital plane is misaligned 
with the stellar rotation axis (e.g. Hebrard et al., 2008; Winn 
ct al. 2009a, 2009b; Pont et al., 2009a, 2009b; Gillon 2009; 
Johnson et al., 2009; Narita et al. 2009). Winn et al. (2009a) 
conclude that a model on misalignment angles which sup- 
poses two planetary populations - one perfectly aligned and 
one isotropically orientated - is consistent with the data. 
They interpret this as being due to two different causes of 
inward planetary migration - one via a gaseous disc and one 
via planet-planet scattering. We have found, however, that 
even gaseous discs can show substantial variation in mis- 
alignment angles depending on the history or accretion and 
more especially on the history of close stellar interactions 
while accretion is still taking place. 

However, it also needs to be borne in mind that the 
accretion of the disc material remaining at the end of our 
calculations can still affect the stellar rotation. A number 
of authors (Cameron & Campbell, 1993; Yi, 1994; Ghosh, 
1995; Armitage & Clarke, 1996) have argued that the rota- 
tion rates of young stars can be regulated by magnetic link- 
age between the star and a surrounding accretion disc. By 
the same argument, the misalignment of spins between disc 
and star could also be regulated by the same mechanism. 
For example the amount of mass AM needed to significantly 
change the stellar spin is of order AAf/M, ~ (kR„ / ' Rm) 2 
where kR* is the star's radius of gyration and Rm is the 
magnetospheric radius. For typical values of k 2 ~ 0.2 (Ar- 
mitage & Clarke, 1996) and Rm/R* ~ 5-10 (Gregory et al., 
2008) we estimate AM/M, ~ 10~ 2 - 10~ 3 . We note further 
that formation of the solar planets requires a 'minimum so- 
lar nebular' mass of around Md/M* > 0.02. Determining 
how much of this mass is actually able to interact with the 
central star and so change its spin will depend on the de- 
tails of the planet formation process, and in particular on 
the timing and extent of the formation of any inner gap 
which would decouple disc and star. 



6 CONCLUSIONS 

We have considered the evolution of the stellar and disc spin 
axes during the formation of a star which is accreting in 
a variable fashion from an inherently chaotic environment. 



To model this process we have modelled the evolution of a 
warped disc which has material added to it in a variable 
fashion. We take the input of mass and angular momentum 
to the disc to that acquired by the 'sink particle' in an SPH 
simulation of the formation of a cluster of stars in from self- 
gravitating turbulent gas. We have noted that many of the 
assumptions we have used cannot be made self-consistent. 
Thus we caution against drawing specific conclusions from 
our analysis. The calculations we show here should be re- 
garded simply as illustrative. Nevertheless there are some 
general points which can be made. 

First, the variability of the direction of the spin of ma- 
terial accreting onto the central protostar implies that the 
internal velocity field of such stars may be more complicated 
than the usual assumption of aligned differential rotation. 

Second the lighter the disc (and the disc can lose mass 
either through rapid accretion or through interaction with a 
nearby protostar) the more able is the accreted material to 
cause the inner disc spin, which we identify with the direc- 
tion of a jet, to vary. 

Third, the final stellar rotation axis and the final disc 
spin axis can be strongly misaligned. However, this occurs 
most strongly when the disc is truncated so that rotation di- 
rection of the final disc material depends simply on what fell 
in last. In this case, although the disc might be misaligned 
with the star, it might not contain enough material to form 
planets. And we have noted that, depending on the details 
of the models, it may be that a misaligned disc which is mas- 
sive enough to form planets, may also be massive enough to 
reduce the misalignment. 

In conclusion, we have shown that it should be pos- 
sible make some deductions about the accretion history of 
a young star by observations of jet direction variability, of 
star/planetary orbit misalignments, and perhaps even of its 
internal rotation structure. It is clear that further work is re- 
quired to improve models of the star formation process with 
better account being taken of the nature of the accretion 
process within the final 100 au or so. 
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APPENDIX: ADDITION AND SHIFTING OF 

MASS 

We discuss here our algorithm for adding material to the 
disc. 

Suppose that the disc grid cells are at radii Ri,i = 1,N 
corresponding to specific angular momenta hi,i = 1,N. In 
a timesetep At we add a mass AM with angular momen- 
tum A J. Thus the specific angular momentum of the added 
material is 

h =4r7, (9) 

AM v ' 

where AJ = |AJ|. Matter is than added to the cells k and 
k + 1 where hk < h < h k +i- 

Then the angular momentum in cell k is incremented 
by an amount 

hk+i — h 



AJ fc = 



-AJ, 



(10) 



hk+i — kk 
and that in cell k + 1 by an amount 

AJ fc+1 = -A— \ AJ. (11) 

tlk+l — Kk 

In this manner angular momentum is exactly conserved. 
However, the amount of mass actually added is now 



AM , = AJk ±1+ AJk^ 
hk+i hk 



(12) 



which needs to be compared with the actual mass to be 
added AM = AJ/h. 

If we write hk+i = h{\ + €2) and hk = h(l — €1), where 



Alignment of disc and stellar rotation axes 



€2 > and ei > then to first order in the small quantities 
e it is simple to show that 

AM' -AM 

AM =™>°- ( 13 ) 

Thus the algorithm for the addition of angular momen- 
tum leads to a systematic over-estimate of the amount of 
mass added of an amount which is on average (i.e. when ei 
and £2 are of comparable magnitude) second order in dR/R. 



